package gbench.sandbox.matlib;

import org.junit.jupiter.api.Test;

import java.text.NumberFormat;
import java.util.*;
import gbench.common.fs.FileSystem;
import gbench.common.matlib.MatlibCanvas.NVec;
import gbench.common.matlib.MatlibCanvas.Range;
import gbench.common.matlib.MatlibCanvas.Seq;
import gbench.common.matlib.MatlibCanvas.XVec;
import gbench.common.matlib.data.DataReader;
import gbench.common.tree.LittleTree.IRecord;

import static gbench.common.matlib.MatlibCanvas.println;
import static gbench.common.matlib.data.Pipe.DUMMY;

/**
 * 广告数据回归分析
 * @author gbench
 *
 */
public class JunitAdvertising {
    
    @Test
    public void foo() {
        final var file = FileSystem.fileOf("advertising.txt", this.getClass()); // 读取广告数据文件
        final var dfm = DataReader.loaddf(file,1);// 读取数据: 0 广告数据, 1 模拟数据, 算法对 广告数据 回归不理想
        final var dfm_train = dfm  // 抽样训练数据
            // .shuffle() // 打乱顺序
            .head(100); // 提取100条数据用作训练
        println("dfm_train 训练数据集合\n", dfm_train);
        
        final var bundle = dfm_train.dfrm(rec -> rec.update((String k, List<String> col) -> { // 列类型转换
            return col.stream().map(IRecord.coerce(Double.class)).toArray();
        })).split(Range.of(4)).fmap(bp -> bp // 裂解的行项变换
            .fmap1(r -> r.arrayOf(Double.class,NVec::new).set(0, 1d)) // 提取一组元素,并在首位设置常量系数1
            .fmap2(r -> r.dbl(0)) // 只提取一个元素,生成一条值向量
        ).unwind(); // 纵向分解为向量束
        final var x_train = bundle._1(); // 训练数据集合
        final var y_train = bundle._2().mutate(NVec::new); // 基准数据集合
        
        // 回归分析
        println("参数thetas的估计值（回归结果）:",regress(x_train,y_train).toString(nf::format));
    }
    
    /**
     * 多元线性回归
     * 
     * @param x_train 训练数据的 集
     * @param y_train 训练数据的 标准值
     */
    @Test
    public NVec regress(final XVec<NVec> x_train,final NVec y_train) {
        
        final var num_variables = x_train.shape().get(1); // 获取训练数据集的宽度
        final var mx = x_train.map(e->e.toArray()).toArray(Double[][]::new);
        final var x_train_t = XVec.OF(Seq.transpose(mx)).fmap(NVec::new); // x_train的转置 也是 XVec<NVec> 类型的结构
        final var thetas = NVec.of(num_variables, i -> 0); // 估计的参数
        final var error_pipe = DUMMY((XVec<NVec>) null) // 输入一条观测行向量组：XVec<NVec> 即源数据矩阵 x_train , 管道类型为: Pipe<XVec<NVec>, NVec>
            .bind(xvec -> xvec.fmap(thetas::dot).mutate(NVec::new)) // x_i_j*theta_j -> y_i
            .bind(yhat -> yhat.sub(y_train)); // 残差计算管道 yhat:[y_i]
        
        println("\nx_train 训练数据x矩阵:\n", x_train.toString(e -> e + "\n"));
        println("\nx_train_t 训练数据x矩阵的转置 :\n", x_train_t.toString(e -> e + "\n")); // 转置是为了方便取列数据 以方便计算：error.dot(x_train_t.get(j)) 即j索引方向的梯度分量
        println("\ny_train 训练数据实际值:\n", y_train.toString(e -> e + "\n"));
        println("参数thetas的估计值（初始状态,默认为0):" + thetas.toString(nf::format));
        
        ///////////////////////////////////////////////////////////////////////////////
        // 训练过程的参数基础设置
        ///////////////////////////////////////////////////////////////////////////////
        
        final var rate = 1e-2; // 学习率-步长
        final var delta = 1e-7; // 小参数,用于防止 分母为0, 即 r.get(j) 为0
        final var eps = 1e-3; // 误差容许范围：最大值的误差累计和
        final var m = x_train.size(); // 数据记录的数量
        final var error = NVec.of(x_train.size(), i -> 0); // 残差向量，可以采用error_pipe.evaluate(x_train);来模拟计算一下error向量长度，但该结果在还会循环体重新计算一次。故此处使用0向量做占位。
        final var r = NVec.of(num_variables, i -> 0); // 梯度积累变量
        
        var i = 0; // 循环变量
        var loss = Double.MAX_VALUE; // loss 是一个error向量函数，把error向量生成一个loss指标。以便做容忍限度eps的检测，即调整thetas使得loss达到容许范围
        // 计算残差 以及 相应的 损失函数 直至 达到 低于 容许范围 eps要求
        while ( (loss = error.assign(error_pipe.evaluate(x_train)).fmapN(e -> e * e).mean()) >= eps ) { // LOSS = MSE（Mean Square Error ）均方误差
            if (i++ % 10000 == 0) println(i, "#", loss, thetas); // 显示残差大小
            // 注意 - rate 是由 error_pipe 的 yhat-y_train来定义的。符号(-)代表与 yhat的增加方向（梯度方向)相反,即 负梯度方向
            thetas.update((j, theta_j) -> {
                final var grad_j = 2 * error.dot(x_train_t.get(j)) / m; // 梯度的 j 分量
                r.idval_add(j, grad_j * grad_j); // 梯度的平方和的累加
                final var delta_theta_j = -rate * grad_j / Math.sqrt(r.get(j) + delta); // AdaGrad算法 计算 更新增量
                return theta_j + delta_theta_j; // 更新参数 theta_j
            }); // thetas 系数重新计算
        } // while
        
        return thetas;
    }

    final NumberFormat nf = NumberFormat.getInstance(); // 数字格式化函数
}
